#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mkl.h>          // Intel MKL头文件
#include <sys/time.h>     // 高精度计时

double get_current_time() {
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + tv.tv_usec * 1e-6;
}

int main() {

    mkl_set_num_threads(1);

    const int m = 512, n = 1024, k = 2048;
    const double flops = 2.0 * m * n * k;  // 总浮点运算次数

    // 使用MKL内存分配（自动64字节对齐）
    double *A = (double*)mkl_malloc(m * n * sizeof(double), 64);
    double *B = (double*)mkl_malloc(n * k * sizeof(double), 64);
    double *C = (double*)mkl_malloc(m * k * sizeof(double), 64);

    // 初始化随机矩阵
    srand(time(NULL));
    for (int i = 0; i < m * n; i++) A[i] = (double)rand() / RAND_MAX;
    for (int i = 0; i < n * k; i++) B[i] = (double)rand() / RAND_MAX;
    for (int i = 0; i < m * k; i++) C[i] = 0.0;

    // 预热运行（消除冷启动影响）
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 
                m, k, n, 1.0, A, n, B, k, 0.0, C, k);

    // 正式计时
    double start = get_current_time();
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 
                m, k, n, 1.0, A, n, B, k, 0.0, C, k);
    double elapsed = get_current_time() - start;

    // 性能计算
    double gflops = (flops / elapsed) / 1e9;

    printf("\n===== MKL矩阵乘法性能 =====\n");
    printf("矩阵维度: %d x %d x %d\n", m, n, k);
    printf("计算时间: %.6f 秒\n", elapsed);
    printf("性能: %.2f GFLOPS\n", gflops);
    printf("使用的函数: cblas_dgemm()\n");
    printf("==========================\n");

    // 释放MKL分配的内存
    mkl_free(A);
    mkl_free(B);
    mkl_free(C);

    return 0;
}